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Abstract 

A  numerically  efficient  mathematical  model  of  a  proton  exchange  membrane  fuel  cell  (PEMFC)  stack  is  presented.  The  aim  of  this  model  is 
to  study  the  dynamic  response  of  a  PEMFC  stack  subjected  to  load  changes  under  the  restriction  of  short  computing  time.  This  restriction  was 
imposed  in  order  for  the  model  to  be  applicable  for  nonlinear  model  predictive  control  (NMPC).  The  dynamic,  non-isothermal  model  is  based  on 
mass  and  energy  balance  equations,  which  are  reduced  to  ordinary  differential  equations  in  time.  The  reduced  equations  are  solved  for  a  single  cell 
and  the  results  are  upscaled  to  describe  the  fuel  cell  stack.  This  approach  makes  our  calculations  computationally  efficient.  We  study  the  feasibility 
of  capturing  water  balance  effects  with  such  a  reduced  model.  Mass  balance  equations  for  water  vapor  and  liquid  water  including  the  phase  change 
as  well  as  a  steady- state  membrane  model  accounting  for  the  electro-osmotic  drag  and  diffusion  of  water  through  the  membrane  are  included. 
Based  on  this  approach  the  model  is  successfully  used  to  predict  critical  operating  conditions  by  monitoring  the  amount  of  liquid  water  in  the  stack 
and  the  stack  impedance.  The  model  and  the  overall  calculation  method  are  validated  using  two  different  load  profiles  on  realistic  time  scales  of  up 
to  30  min.  The  simulation  results  are  used  to  clarify  the  measured  characteristics  of  the  stack  temperature  and  the  stack  voltage,  which  has  rarely 
been  done  on  such  long  time  scales.  In  addition,  a  discussion  of  the  influence  of  flooding  and  dry-out  on  the  stack  voltage  is  included.  The  modeling 
approach  proves  to  be  computationally  efficient:  an  operating  time  of  0.5  h  is  simulated  in  less  than  1  s,  while  still  showing  sufficient  accuracy. 

©  2008  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

In  recent  years  the  interest  in  using  hydrogen  fuel  cells  as 
power  supply  for  portable  electronics  has  grown  substantially. 
Compared  to  batteries  fuel  cell  systems  can  provide  a  higher 
energy  density  and  instantaneous  refilling  while  avoiding  the 
problem  of  self-discharge.  However,  the  use  of  fuel  cells  as 
power  supply  for  electronic  products  is  challenging  because  the 
power  demand  of  these  applications  fluctuates.  Due  to  the  lim¬ 
ited  space  in  portable  electronics  the  stack  can  in  many  cases 
not  be  buffered  by  a  battery.  Thus,  the  fuel  cell  does  not  usually 
operate  at  steady- state.  A  solid  understanding  of  the  dynamic 
response  of  a  proton  exchange  membrane  fuel  cell  (PEMFC) 
under  load  changes  is  crucial  for  reliable  and  optimized  opera¬ 
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tion  [1].  The  dynamic  behavior  of  a  fuel  cell  is  a  highly  complex 
phenomenon,  as  it  involves  different  length  and  time  scales.  The 
power  of  a  PEMFC  also  strongly  depends  on  operating  condi¬ 
tions  such  as  flow  rates,  relative  humidity  and  temperature  of 
the  gases  as  well  as  ambient  temperature.  Mathematical  mod¬ 
eling  is  a  powerful  tool  for  understanding  and  handling  this 
complexity.  Nonlinear  model  predictive  control  (NMPC)  and 
online  optimization  of  dynamic  processes  have  attracted  increas¬ 
ing  attention  over  the  past  decade,  see  e.g.  Ref.  [2].  In  contrast 
to  empirical  control  strategies  based  on  experimental  observa¬ 
tions  and  extensive  testing,  a  model-based  control  allows  faster 
system  development  and  optimal  system  operation  over  a  wide 
range  of  operating  conditions.  As  a  prerequisite,  NMPC  requires 
detailed  nonlinear  process  models. 

A  considerable  amount  of  work  has  been  done  thus  far  to 
model  PEMFCs  [3,4].  Most  of  the  models  are  steady-state, 
see  for  example  Refs.  [5-7].  Less  work  has  been  published 
on  dynamic  fuel  cell  modeling.  Amphlett  et  al.  [8]  modeled 
the  behavior  of  the  stack  temperature  and  the  voltage  during 
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start-up,  shut-down  and  load- step.  In  their  model  only  the  energy 
balance  of  the  solid  is  modeled  dynamically  whereas  all  other 
equations  are  assumed  to  be  at  quasi-steady  state  for  a  given 
solid  temperature.  Lee  et  al.  [9]  used  an  object-oriented  approach 
based  on  stationary  equations.  Dynamic  profiles  are  created  by 
calculating  the  quasi- stationary  solution  variables  for  each  time- 
step.  Ceraolo  et  al.  [10]  used  an  isothermal  model  to  simulate 
the  dynamic  behavior  of  the  cell  voltage  to  a  load  change  on  a 
time-scale  of  seconds.  Their  model  was  extended  to  account  for 
non-isothermal  conditions  by  Shan  et  al.  [11].  In  recent  years 
several  authors  have  presented  dynamic  models  using  a  similar 
approach  as  Amphlett  et  al.  [8].  Golbert  et  al.  [12]  developed  a 
transient  along-the-channel  model  for  control  purposes,  which 
includes  mass  balances  of  liquid  water  and  water  vapor.  Yu  et 
al.  [13]  presented  another  extension  of  the  model  of  Amphlett, 
which  accounts  for  the  influence  of  latent  heat  on  the  energy 
balance.  Pathapati  et  al.  [14]  included  dynamic  mass  balance 
equations  and  energy  balance  equations  for  the  gases.  In  addi¬ 
tion,  a  term  to  account  for  the  double  layer  capacity  is  presented. 
The  influence  of  flooding  on  the  dynamic  behavior  of  the  stack 
voltage  under  isothermal  conditions  was  modeled  by  McKay  et 
al.  [15].  All  of  these  models  are  reduced  in  terms  of  dimension¬ 
ality  and  comprehensiveness.  In  recent  years  several  research 
groups  have  published  valuable  in-depth  analyses  on  the  tran¬ 
sient  behavior  of  PEM  fuel  cells  using  commercial  software  tools 
either  based  on  computational  fluid  dynamics,  e.g.  Refs.  [16-21] 
or  on  finite-element,  multiphysics  simulation  approaches,  e.g. 
Refs.  [22-25].  These  studies  help  in  understanding  the  funda¬ 
mental  physical  processes  and  interactions  within  the  fuel  cell. 
However,  due  to  the  massive  computational  effort  required  they 
are  not  suitable  for  online  control. 

This  work  presents  a  dynamic  model  approach  for  portable 
fuel  cell  stacks.  The  aim  of  our  work  is  to  model  the  dynamic 
behavior  of  a  portable  PEM  fuel  cell  stack  on  relevant  time 
scales  for  technical  applications  under  the  restriction  of  keeping 
the  computing  time  short.  Therefore,  a  reasonable  compromise 
between  physical  accuracy  and  numerical  efficiency  is  found 
which  makes  the  model  suitable  for  NMPC.  Despite  of  the 
number  of  very  valuable  contributions  to  the  dynamic  fuel  cell 
modeling,  a  validated  model  approach,  which  meets  these  needs 
was  not  found  in  the  literature.  Many  of  the  existing  models 
require  massive  computational  effort  either  in  terms  of  memory 
usage  or  computing  time  or  both,  e.g.  Refs.  [16-21,23-25].  Most 
of  the  reduced  models  in  the  literature  are  either  not  designed  for 
NMPC  purposes,  e.g.  Refs.  [9,10]  or  are  not  validated  against 
experimental  data  of  a  PEMFC  stack  on  realistic  time  scales, 
e.g.  Refs.  [11-13].  The  reduced  model  presented  here  is  vali¬ 
dated  against  experimental  data  of  a  PEMFC  stack  for  different 
load  profiles  on  realistic  time  scales  of  up  to  30  min.  The  model 
validation  study  does  also  include  an  analysis  of  the  character¬ 
istics  of  the  stack  in  critical  states  of  operation.  The  model  is 
non-isothermal  and  considers  the  mass  transfer  and  the  electro¬ 
chemical  reactions.  Moreover,  the  model  accounts  for  both  water 
vapor  and  liquid  water  and  for  the  phase  transition.  In  contrast  to 
previous  modeling  studies  the  average  liquid  water  concentra¬ 
tion  is  used  to  predict  flooding  of  the  fuel  cell  based  on  the  liquid 
water  balance.  This  can  be  useful  for  improved  NMPC  algo¬ 


rithms.  To  ease  the  transfer  to  different  stacks  we  describe  the 
methods  of  parameter  identification  in  detail.  In  order  to  meet  the 
challenge  of  realizing  sufficiently  exact  modeling  results  with 
short  calculation  time  some  strong  simplifications  are  made, 
which  are  justified  by  the  good  agreement  between  simulation 
and  experimental  results. 

2.  Model  description 

2.7.  Modeling  approach 

The  PEM  fuel  cell  stack  model  presented  here  is  dynamic  and 
non-isothermal.  The  model  is  based  on  transient  energy  and  mass 
balance  equations,  a  membrane  model  and  an  electrical  model 
based  on  the  tafel  equation.  Convective  heat  and  mass  transfer 
within  the  stack  are  accounted  for  dynamically.  A  mass  balance 
of  water  in  the  liquid  and  vapor  phase  is  included.  Condensation 
and  evaporation  in  the  channels  as  well  as  water  generation  at 
the  cathode  are  accounted  for.  The  membrane  model  is  steady 
state  and  accounts  for  the  electro-osmotic  drag  and  back  dif¬ 
fusion  of  water.  The  steady-state  electrical  model  incorporates 
the  influence  of  pressure  and  temperature  changes  as  well  as  the 
voltage  drop  due  to  activation  and  ohmic  losses. 

Fuel  cell  stacks  are  commonly  characterized  by  measuring 
the  time  evolution  of  the  stack  voltage  and  the  stack  tempera¬ 
ture  subject  to  specific  operating  conditions  like  power  demand, 
ambient  and  gas  temperature  and  gas  humidity.  In  integrated  fuel 
cell  systems  the  stack  voltage,  the  stack  temperature  and  the  gas 
flow  rates  are  usually  monitored.  The  stack  model  presented  here 
allows  the  simulation  of  the  most  important  parameters  for  the 
operation  of  a  PEMFC  stack.  The  model  considers  four  mon¬ 
itoring  points  for  the  mass  and  heat  balance  of  the  gases:  the 
inlet  and  outlet  of  the  stack  on  the  anode  and  cathode  side.  At 
these  points,  gas  temperatures  and  molar  fluxes  of  the  different 
species  are  considered.  The  operating  conditions  at  the  stack 
inlets  are  used  as  input  values  of  the  model.  Furthermore,  the 
current  density  is  an  input  variable.  Based  on  the  values  of  the 
operating  conditions,  the  model  predicts  the  molar  fluxes  and  gas 
temperatures  at  the  stack  outlets,  the  average  temperature  of  the 
solid  material,  the  stack  voltage  and  the  average  concentration 
of  liquid  water  in  the  stack.  Fig.  1  shows  the  solution  variables 
and  operating  conditions  and  illustrates  that  these  parameters 
are  accessible  even  in  a  fully  integrated  fuel  cell  system.  For  the 
stack  temperature,  four  monitoring  points  are  indicated  from 
which  the  average  temperature  is  calculated. 

The  derivation  of  the  model  equations  is  done  in  a  bottom-up 
approach  that  can  be  split  up  into  three  steps.  First,  we  con¬ 
sider  one  particular  cell  of  the  stack.  Balance  equations  are  set 
up  for  a  representative  elementary  volume  of  the  cell  (REVi) 
and  for  a  representative  elementary  channel  (REC)  in  the  cell 
(REV2).  The  model  is  reduced  to  one  geometrical  dimension,  the 
direction  along  the  channel  length.  Second,  the  time-dependent 
balance  equations  for  the  REVs  are  integrated  along  the  channel 
length.  For  a  complete  description  of  one  cell  all  gas  channels 
are  assumed  to  behave  like  the  REC.  Third,  the  stack  is  modeled 
as  several  coupled  cell  modules.  The  approach  and  the  model 
geometry  is  illustrated  in  Figs.  1-3.  The  model  assumptions  are 
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Fig.  1.  3D-View  of  a  PEMFC  stack  consisting  of  six  cells.  The 
membrane-electrode  assemblies  between  the  bipolar  plates  are  indicated  in  dark 
gray.  The  measurement  points  that  were  used  for  the  characterization  of  the  stack 
with  respect  to  thermal  and  mass  balance  are  indicated.  The  values  obtained  at 
these  measurement  points  correspond  to  the  input  and  solution  variables  of  the 
model.  The  cut  for  Fig.  2  is  indicated  in  light  gray. 

motivated  by  the  preservation  of  computational  efficiency  which 
is  required  for  the  simulation  of  realistic  load  cycles  in  system 
simulation  and  for  the  application  in  NMPC.  In  contrast,  more 
detailed  model  approaches,  e.g.  Refs.  [26,27]  require  massive 
computational  effort  either  in  terms  of  memory  usage  or  com¬ 
puting  time  or  both.  The  validity  of  our  approach  is  shown  by  a 
comparison  of  simulation  and  experimental  results  for  different 
realistic  load  profiles. 

2.2.  Model  assumptions 

The  model  assumptions  are  listed  below.  The  axis  system  is 
illustrated  in  Fig.  2. 

•  The  gas  channels  are  treated  as  plug-flow  reactors.  This 
approach  is  used  by  several  other  authors,  e.g.  Refs. 
[6,10,12,28-30]  and  is  generally  accepted.  A  good  overview 
of  approaches  used  for  the  modeling  of  gas  flow  in  the  chan¬ 
nels  is  given  in  Refs.  [3,4].  A  uniform  velocity  profile  and 


Fig.  2.  3D-View  of  the  cell  model  corresponding  to  the  cut  in  Fig.  1,  indicating 
the  cross-sectional  area  of  the  cell  As  and  of  the  gas  channels  Aa  and  Ac,  the 
representative  elementary  volume  REVi  as  well  as  the  cut  for  Fig.  3. 


Fig.  3.  Sideview  of  a  cell  along  the  cut  indicated  in  Fig.  2,  illustrating  the  direc¬ 
tion  of  the  molar  flow  Naj  of  species  i  in  the  anode  gas  channel  and  the  molar 
flow  Ncj  of  species  j  in  the  cathode  gas  channel  from  the  channel  inlets  to  the 
channel  outlets.  The  representative  elementary  volume  REV2  is  exemplified  for 
the  anode  gas  channel. 

a  complete  and  instantaneous  mixing  perpendicular  to  the 
direction  of  the  gas  flow  is  assumed.  Thus,  the  gas  velocity 
and  concentrations  within  a  gas  channel  are  uniform  in  the  y- 
and  z-direction.  The  underlying  assumption  of  laminar  flow 
is  widely  used  in  fuel  cell  modeling,  e.g.  Refs.  [16,21,31-33], 
and  was  also  satisfied  by  a  Reynolds  number  calculation. 

•  The  stack  temperature  is  constant  in  the  y-  and  z-direction. 

•  Liquid  water  is  assumed  to  exist  in  the  form  of  small  droplets 
on  the  surface  of  the  gas  channels  [5].  As  heat  exchange  with 
the  solid  material  is  much  faster  than  with  the  gas,  liquid  water 
is  further  assumed  to  have  stack  temperature. 

•  The  volume  of  the  liquid  water  is  assumed  negligible.  Hence, 
it  has  no  influence  on  the  gas  transport  in  the  channels. 

•  Cross-over  of  nitrogen,  hydrogen  and  oxygen  through  the 
membrane  is  neglected. 

•  The  electrochemical  reactions  and  the  transport  processes  are 
assumed  to  be  homogeneous  throughout  the  stack. 

•  The  mass  transport  resistance  of  the  gas  diffusion  layer  is 
neglected. 

•  The  electrical  and  thermal  contact  resistances  are  neglected. 

Computational  efficiency  is  mainly  achieved  by  the  following 
model  characteristics:  in  order  to  reduce  the  number  of  mod¬ 
ules  to  be  calculated  only  one  cell  is  modeled.  The  results  are 
upscaled  according  to  the  number  of  cells  in  the  stack.  As  an 
adequate  modeling  of  the  highly  complex  phenomena  in  the 
GDL  and  the  membrane  would  require  massive  computational 
effort,  the  GDL  is  neglected  and  the  membrane  is  assumed  to 
be  at  steady-state.  The  spatially  resolved  characteristics  of  the 
fuel  cell  stack  are  neglected  through  the  reduction  of  the  model 
equations  to  ordinary  differential  equations  in  time. 

3.  Mathematical  model 

In  the  first  step  of  the  model  derivation  equations  for  energy 
and  mass  balance  are  formulated  for  a  representative  elementary 
volume  of  the  solid  material  (REV  1 )  and  of  a  representative  gas 
channel  (REV2).  As  illustrated  in  Figs.  2  and  3,  REVi  corre¬ 
sponds  to  a  section  of  a  representative  cell,  whereas  REV2  is 
a  channel  section  of  a  representative  elementary  gas  channel. 
Mass  balance  equations  are  set  up  for  hydrogen  H2,  oxygen  O2, 
nitrogen  N2,  water  vapor  H20v,  and  liquid  water  H201in  REV2. 
The  structure  of  the  mass  balance  equations  is  identical  for  all 
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species: 

dfCk,i  =  —(V  •  Mfrj)  +  Rkj,  (1) 

where  the  change  in  concentration  ck,i  of  species  i  in  volume 
element  k  is  given  by  the  divergence  of  the  molar  flow  rate  Mkj 
and  by  appropriate  rates  of  production  and  consumption  Rkj. 
Energy  balance  equations  are  formulated  for  the  gas  channels  of 
anode  and  cathode  as  well  as  for  the  solid  material.  Their  basic 
structure  is 

dtUk  =  -(V  •  UkvKi)  -  V  •  (jcVT)  +  Sk,  (2) 

where  the  change  of  internal  energy  Uk  in  volume  element  k  is 
given  by  convective  transport  of  energy  associated  with  the  mass 
average  velocity  vkj,  transfer  by  heat  conduction  with  the  heat 
conduction  coefficient  tc,  and  energy  production  or  consumption 
with  the  rate  Sk. 

As  mentioned  previously,  we  assume  plug-flow  conditions 
for  the  gas  channels.  Hence,  Eq.  (1)  can  easily  be  reduced  to  the 
following  equation: 

dt£k,i(.x,  t)Ak  —  dxNkj(x,  t )  +  Rk,i(x,  t)Ak,  (3) 


where  Ak  is  the  cross-sectional  area  of  REV2  and  Nkj  is  the 
molar  flux  of  species  i  along  the  REV.  Analogously,  Eq.  (2)  for 
the  energy  balance  of  both  REVs  is  reduced  to 

dtuk(x,  t)=  -  dx(Uk(x,  t)vkJ(x,  t))-Kd2xT(x,  t)+Sk(x,  t).  (4) 

As  a  consequence  of  this  reduction,  the  mathematical  model 
contains  several  effective  parameters  that  include  spatially  dis¬ 
tributed  effects  in  the  y-z- -plane.  In  Section  6  it  is  described  how 
these  parameters  are  obtained.  Below,  the  different  balance  equa¬ 
tions  are  presented  in  the  concise  form  corresponding  to  Eqs.  (3) 
and  (4),  which  is  the  form  after  step  one  of  the  derivation.  Steps 
two  and  three  are  illustrated  in  Section  4.  The  parameters  and 
symbols  used  in  the  model  equations  are  listed  in  Tables  1  and  2. 
Table  3  contains  a  list  of  subscripts  and  superscripts. 

3.1.  Energy  balance 

3.1.1.  Energy  balance,  gas  channels 

The  following  continuity  equation  describes  the  energy  bal¬ 
ance  of  the  gaseous  species  in  REV2.  The  considered  species  are 
i  =  H2,  H2Ov  for  the  anode  and  i  =  O2,  N2,  H20v  for  the  cath¬ 
ode.  Liquid  water  is  assumed  to  exist  in  form  of  small  droplets 
at  the  surface  of  the  channels.  The  liquid  water  is  assumed  to  be 


Table  1 

List  of  parameters  and  constants 


Symbol 

Explanation 

Value 

Reference 

Aa 

Cross-sectional  area  of  anode  gas  channel 

(9  ±2)  x  10-7  m2 

Meas. 

A-a , sg 

Heat  exchange  area  per  unit  length  between  solid  and  gas,  anode  side 

(4.2  ±  0.4)  x  10“3  m 

Meas. 

Ac 

Cross-sectional  area  of  cathode  gas  channel 

(12  ±2)  x  10“7m2 

Meas. 

Ac,  sg 

Heat  exchange  area  per  unit  length  between  solid  and  gas,  cathode  side 

(4.6  ±0.4)  x  10“3m 

Meas. 

As 

Cross-sectional  area  of  cell 

(5.7  ±0.9)  x  10-5  m2 

Meas. 

Ass 

Heat  exchange  area  between  solid  and  surroundings  per  unit  length 

(1.9  ±0.2)  x  10-2  m 

Meas. 

Ch2 

Heat  capacity  of  H2 

28.8  J  mol-1  K-1 

[37] 

ch2o1 

Heat  capacity  of  H2O1 

ys.sjmor'K-1 

[37] 

ch2ov 

Heat  capacity  of  H20v 

SS^Jraor'K-1 

[37] 

Cn2 

Heat  capacity  of  N2 

29.1  J  mol- ‘K"1 

[37] 

Co2 

Heat  capacity  of  O2 

29.3  Jmol-1  K-1 

[37] 

Cs 

Specific  heat  capacity  of  stack  solid  material 

(770  ±  40)  J  kg"1  K-> 

Meas. 

Dref 

m,H20 

Diffusion  coefficient  of  water  in  membrane  for  standard  conditions 

5.5  x  10-um2s-1 

[36] 

dy 

Scaling  factor  for  the  channel  width 

2.2  ±0.2 

Meas. 

F 

Faraday  constant 

96,  485  C  mol-1 

[38] 

K 

Height  of  anode  gas  channel 

(6  ±  1)  x  10-4m 

Meas. 

hc 

Height  of  cathode  gas  channel 

(8  ±  1)  x  10-4  m 

Meas. 

A  H\ap 

Enthalpy  of  phase  transition  of  water 

44,  000 Jmol-1 

[37] 

^phase 

Condensation  rate  constant 

100  s"1 

[36] 

^m,p 

Permeability  of  water  in  membrane 

1.58  x  10-18  m2 

[36] 

Le  ff 

Effective  channel  length 

(465  ±  2)  x  10“3  m 

Meas. 

n  cells 

Number  of  cells 

6 

Meas. 

^chan 

Number  of  channels  in  each  flow  field 

2 

Meas. 

R 

Ideal  gas  constant 

8.314JK-1  mol-1 

[38] 

ASa 

Molar  entropy  of  reaction  at  anode 

0.104  JKT1  mol-1 

[39] 

ASC 

Molar  entropy  of  reaction  at  cathode 

— 326.36JK-1  mol-1 

[39] 

2^0 

Reference  temperature 

298  K 

Meas. 

fin 

Thickness  of  membrane 

(250  ±  1)  x  10“7  m 

Meas. 

U° 

*^oc 

Open-circuit  voltage  for  standard  conditions 

1.23  V 

Meas. 

Esg 

Heat  transfer  coefficient  between  solid  material  and  gas 

25  W  m-2  K-1 

[5] 

Ess 

Heat  transfer  coefficient  between  solid  material  and  surroundings 

4.7  ±  0.9  Wm- 2  K"1 

Meas. 

Wq 

Width  of  gas  channel 

(1.5  ±0.1)  x  10“3m 

Meas. 

Pn2ol 

Viscosity  of  H2O1 

3.56  x  10-4Pas 

[37] 

Ps 

Average  density  of  stack  material 

2500  ±  100  kg  m"3 

Meas. 
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Table  2 


List  of  symbols 


Symbol 

Explanation 

Unit 

A  k 

Cross-sectional  area  of  REV2 

m2 

Cross-sectional  area  of  gas  channel  in  electrode  q 

m2 

«H2Ov 

Activity  of  water  vapor 

- 

Ck 

Constant  with  index  k 

- 

Ck,i 

Concentration  of  species  i  in  volume  element  k 

mol  m-3 

Cq,i 

Concentration  of  species  i  in  gas  channel  q 

mol  m-3 

An,H20 

Diffusion  coefficient  of  H2O  in  the  membrane 

m2s-1 

I 

Current 

A 

j 

Current  density 

Am-2 

70, c 

Cathodic  exchange  current  density 

Am-2 

Mk,i 

Molar  flow  rate  of  species  i  in  volume  element  k 

mols-1  m-2 

Nq,i 

Molar  flow  of  species  i  in  gas  channel  q 

mols-1 

Nk,i 

Molar  flow  of  species  i  in  REV2 

mol  s_1 

.  ,in/out,rec 

iyq3 

Molar  flow  of  species  i  into/out  of  the  REC  of 
electrode  q 

mol  s_1 

Arin,cell 

1Sq,i 

Molar  flow  of  species  i  into  the  gas  inlet  of  the 
cell  on  the  q  side 

mols-1 

at  out, cel  1 
^  q,i 

Molar  flow  of  species  i  out  of  the  gas  outlet  of  the 
cell  on  the  q  side 

mols-1 

Arin, stack 
^  q,i 

Molar  flow  of  species  i  into  the  gas  inlet  of  the 
stack  on  the  q  side 

mols-1 

Arout, stack 

nq,i 

Molar  flow  of  species  i  out  of  the  gas  outlet  of  the 
stack  on  the  q  side 

mol  s_1 

At q,  phase 

Rate  of  phase  change  in  electrode  q 

mol  s_1  m_l 

^drag 

Electro-osmotic  drag  coefficient 

- 

p. 

Average  pressure  in  gas  channel  q 

Pa 

Pq,  sat 

Saturation  pressure  in  gas  channel  q 

Pa 

Pq,i 

Partial  pressure  of  species  i  in  electrode  q 

Pa 

P? 

Partial  pressure  of  species  i  at  reference 
conditions 

Pa 

Rk,i 

Rate  of  consumption  or  production  of  species  i  in 
volume  element  k 

mol  m-3  s_1 

Sk 

Rate  of  energy  conversion  in  volume  element  k 

W  m-3 

T 

Temperature 

K 

Ta,m 

Gas  temperature  at  anode  gas  inlet 

K 

Tl,OUt 

Gas  temperature  at  anode  gas  outlet 

K 

T;,in 

Gas  temperature  at  cathode  gas  inlet 

K 

T,out 

Gas  temperature  at  cathode  gas  outlet 

K 

y^rec 
q,  in/out 

Gas  temperature  at  gas  outlet/inlet  of  REC 

K 

Tq 

Gas  temperature  in  gas  channel  q 

K 

Ts 

Stack  temperature 

K 

^s,av 

Average  stack  temperature 

K 

Tur 

Temperature  of  surroundings 

K 

t 

Time 

s 

t^cell 

Cell  voltage 

V 

Uk,i 

Internal  energy  of  species  i  in  volume  element  k 

Jm-3 

Uoc 

Open-circuit  voltage 

V 

U stack 

Stack  voltage 

V 

Vk,i 

Velocity  of  species  i  in  volume  element  k 

ms-1 

Vq,i 

Velocity  of  species  i  in  gas  channel  q 

ms-1 

Z 

Number  of  electrons  exchanged  in  reaction 

- 

a 

Symmetry  factor 

- 

a  net 

Net  water  migration  coefficient 

- 

Id 

Activation  losses  in  cathodic  reaction 

V 

hoc 

Losses  in  open-circuit  voltage 

V 

hohm 

Ohmic  losses 

V 

Ks 

Heat  conduction  coefficient  of  solid  material 

W  m_1  K_1 

km 

Membrane  water  content 

- 

V/ 

Stoichiometry  factor  of  species  i 

- 

Membrane  conductivity 

Sm"1 

Table  3 


Subscripts  and  superscripts 


Symbol 

Explanation 

0 

Reference  conditions 

a 

Anode  side 

act 

Active  area 

av 

Average 

c 

Cathode  side 

cell 

One  cell  in  the  stack 

chan 

Channel 

drag 

Electro-osmotic  drag 

eff 

Effective 

H2O1 

Liquid  water 

H2Ov 

Water  vapor 

i 

Species  index 

in 

Gas  inlet 

J 

Species  index 

k 

Volume  element 

1 

Liquid  phase 

m 

Membrane 

oc 

Open-circuit 

out 

Gas  outlet 

phase 

Phase  transition 

9 

Anode  or  cathode 

rec 

Representative  elementary  channel 

ref 

Reference 

s,  stack 

Of  the  stack 

sg 

Exchange  between  solid  and  gas 

ss 

Exchange  between  solid  and  surroundings 

sur 

Surroundings 

V 

Vapor  phase 

at  stack  temperature.  The  energy  balance  of  the  gas  reads 

^  ^t[Cgj(x,  t)Ci  Tq(x,  t)Aq] 
i 

=  Tgj(x,  t)CiTq(x,  AqSgUSg[Ts(x,  t)—Tq(x,  t)\ 

i 

—  A//vapA^^r?phase(-^!>  0>  (5) 

where  q  =  a,  c  denotes  the  anode  and  the  cathode  side,  respec¬ 
tively.  i  is  the  species  index,  cqj  is  the  concentration  of  species 
i  in  the  channel  q ,  Q  is  the  heat  capacity,  Tq  is  the  gas  tem¬ 
perature  and  Aq  is  the  cross-sectional  area  of  the  REC,  which 
equals  A&.  Nqj  is  the  molar  flux  of  species  i  along  the  channel 
q.  Aq  Sg  is  the  heat  exchange  area  between  gas  and  solid  mate¬ 
rial,  Us g  is  the  corresponding  heat  exchange  coefficient,  Ts  is 
the  stack  temperature.  A//vap  is  the  enthalpy  of  evaporation  and 
Nq, phase  is  the  rate  of  the  phase  change  per  unit  length.  The  term 
on  the  left-hand  side  of  the  equation  describes  the  change  of 
internal  energy  in  a  volume  element  of  unit  length.  The  terms 
on  the  right-hand  side  describe  from  left  to  right:  the  transfer  of 
internal  energy  by  convection;  the  heat  transfer  between  gas  and 
solid  stack  material  and  the  heat  consumption  or  production  due 
to  evaporation  and  condensation  of  water. 

3.1.2.  Energy  balance,  solid  material 

The  energy  balance  for  the  solid  material  in  REVi  links 
the  processes  in  the  gas  channels  of  the  anode  and  cathode. 
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It  reads 


equation  of  water  vapor  in  REV2  on  the  cathode  side  reads 


psCsAM(x,  t ) 

—  +^chan  A /Aap[  phased  0+^c, phased  01 

[/A  Sa  A  Sc\  c 

—nch<mdyWq  I  J  0-^oc-hDl“^ohm  j(x,  t) 

— ^chan^sgt^a,  sg(Ts(x,  t )  —  Ta(x,  t ))  +  Ac?  Sg(7's(x,  t ) 

~Tc(x,  0)]  +  Assf/ss[Tsur(x)  —  Ts(x,  t)\  +  AsKsdxTs(x,  t ) 

— ^chan9x([^Va  jj2qi(x,  ^)  +  jj2qi(x,  £)]Cjj2qi  Ts(x,  £))>  (6) 


where  ps  is  the  density  of  the  solid  material,  Cs  is  its  heat  capac¬ 
ity,  As  is  the  cross-sectional  area  of  the  cell  and  nc han  is  the 
number  of  channels  in  each  flow-field.  dy  is  a  scaling  coeffi¬ 
cient,  which  takes  into  account  the  enlargement  of  the  contact 
area  between  gas  and  membrane  due  to  diffusion  under  land 
areas  of  the  flow-field.  wq  is  the  channel  width,  A  Sq  is  the 
entropy  of  reaction  in  electrode  q  and  F  is  the  Faraday  con¬ 
stant.  The  open-circuit  losses  are  r]0 c,  the  activation  losses  in 
the  cathode  are  rfD  and  the  ohmic  losses  are  r]ohm-j  denotes  the 
current  density,  Ass  is  the  heat  exchange  area  between  stack  and 
surroundings,  Uss  is  the  heat  transfer  coefficient,  Tsur  denotes 
the  temperature  of  the  surroundings  and  ks  is  the  heat  conduc¬ 
tion  coefficient.  The  change  of  internal  energy  in  the  solid  stack 
material  in  REV  1 ,  which  is  described  by  the  term  on  the  left-hand 
side,  is  given  by  the  following  source  and  sink  terms  on  the  right- 
hand  side  of  the  equation:  (a)  heat  generation  and  consumption 
due  to  condensation  or  evaporation  of  water  in  anode  or  cath¬ 
ode;  (b)  heat  generation  due  to  activation  energy  and  irreversible 
losses;  (c)  heat  transfer  between  bulk  material  and  gases  of 
anode  and  cathode  side;  (d)  heat  transfer  between  solid  material 
and  surroundings;  (e)  heat  conduction  driven  by  a  temperature 
gradient  within  the  stack;  (f)  heat  transferred  by  convection 
of  liquid  water.  As  a  cell  contains  several  channels  the  terms 
(a)-(c)  and  (f)  need  to  be  multiplied  with  the  number  of  channels 

^chan- 


9/[cc,H2Ov(-*,  *Mc]  =  —  9*[Nc,H2OvO>  0]  -  Nc,phase(*>  0 
dyU)qj(x,  t)  dywqanetj(x,  t) 

+  2F  + - F - ’  (8) 

where  cc?h2ov  is  the  concentration  and  AfCH2ov  is  the  molar  flow 
of  water  vapor  in  REV2.  The  net  water  migration  coefficient  anet 
describes  the  net-number  of  water  molecules  carried  through  the 
membrane  per  proton  (Eq.  (15)).  The  change  in  concentration 
of  water  vapor  is  balanced  by:  the  convective  transport  of  water 
vapor;  the  condensation  or  evaporation  of  water;  the  generation 
of  water  in  the  electrochemical  reaction  and  the  transport  of 
water  through  the  membrane.  The  mass  balance  equation  for 
liquid  water  in  REV2  on  the  cathode  side  is 

0-^c]  =  —  ^X[AC  jj2qi(x,  ^  A^phaseU,  0»  (9) 

where  cQ  H20i  is  the  concentration  and  Nc  H20i  is  the  molar  flow 
of  liquid  water.  The  concentration  of  liquid  water  changes  due 
to  convective  transport  of  liquid  water  along  the  channel  and 
condensation  or  evaporation  of  water.  The  amount  of  water 
that  condenses  or  evaporates  in  a  volume  element  is  modeled 
according  to  Golbert  and  Lewin  [12]: 

V:, phase (x,  0=^7^(Pc, H20»(*.  t)~  Pc^t(Tc(x,  t))),  (10) 
RTc(x,  t) 

where  &phase  is  the  condensation  rate  constant,  hc  is  the  channel 
height,  R  is  the  ideal  gas  constant,  pc,u2ow  denotes  the  partial 
pressure  of  water  vapor  in  the  cathode  gas  channel  and  Pc,sat  is 
the  saturation  pressure  of  water  vapor,  respectively.  It  is  assumed 
that  water  condenses  when  pc,u2ow  >  ^c,sat  and  that  existing 
liquid  water  evaporates  when  pc,u2ow  <  ^c,sat- 


3.3.  Electrical  model 


3.2.  Mass  balance 

For  the  mass  balances  of  REV2  the  mass  transport  along  the 
channel  by  convection,  the  fuel  consumption  of  H2  and  O2,  the 
production  of  H20v  in  the  cathode  side  reaction,  the  evapora¬ 
tion  and  condensation  of  water  as  well  as  the  transport  of  water 
vapor  through  the  membrane  are  taken  into  account.  Below  only 
the  mass  balance  equations  for  the  species  on  the  cathode  side 
are  given.  The  corresponding  equations  for  the  anode  follow  by 
analogy.  The  mass  balance  equation  of  oxygen  in  REV2  on  the 
cathode  side  is 

dvwaI(x,  t ) 

d,[Cc, 02(X,  t)Ac ]  =  -dx[Nc, o2(x,  t)]  -  y-  ^  (7) 

where  cc?o2  is  the  concentration  and  Afc?o2  is  the  molar  flow  of 
oxygen  in  REV2.  The  concentration  of  oxygen  changes  due  to 
convective  transport  along  the  channel  and  due  to  consumption 
of  oxygen  in  the  electrochemical  reaction.  The  mass  balance 


The  cell  potential  C/ce  11  is  calculated  from  the  following  equa¬ 
tion  [12]: 

^cell  —  bl0c  ~  0oc  —  \0d\  —  ^ohm 

tt  1  ( j(f)h n  A 

=  Uoc-rioc - =ln  — - ,  ^  (11) 

uzF  V  jo,c  )  crm(x,  T,  A ,m(0) 

where  Uoc  denotes  the  open-circuit  voltage,  which  is  calculated 
by  the  Nernst-equation  given  in  Appendix  A.  a  is  the  symme¬ 
try  factor  and  z  is  the  number  of  electrons  exchanged  in  the 
reaction.  jo;C  is  the  cathodic  exchange  current  density,  tm  is  the 
thickness  of  the  membrane  and  crm  is  the  membrane  conductiv¬ 
ity.  The  membrane  conductivity  crm  depends  on  the  membrane 
water  content  A.m.  To  account  for  this,  is  calculated  for  each 
time-step.  am  is  calculated  according  to  Springer  et  al.  [34]: 

orm=(0.00514A.m— 0.00326) exp  fl268  •  (12) 
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3.4.  Membrane  model 


Integration  under  use  of  the  mean  value  theorem  leads  to 


Am  — 


The  membrane  water  content  Am  is  defined  as 


lm,S03 


(13) 


where  nm,H2o  is  the  number  of  water  molecules  and  ftm,so3  is 
the  number  of  SO3 -groups  in  the  membrane.  The  dependency 
of  Am  on  the  water  vapor  activity  ^h2ov  is  modeled  according  to 
[34]: 


Am  — 


The  net  water  migration  coefficient  is  given  by  [34] : 

F  cc,H2Ov  _  ca,H2Ov 

C^net  —  drag  .  H20 


( (0.043+17.8aH2ov— 39.85a^20,+36ag20,) 

for  <2h2ov<1 

-l  r 

\(14+1.4(aH2ov-l)) 

for  ftH2ov>l- 

II 

Tfi^ 

_  i 

J  fin 

cc,H2Ov  +  ca,H2Ov  Am,p  F  Pc,U2Ow  ~  /7a,H2Ov 

fin 


(15) 

^h2o  lJ 

where  ft  drag  is  the  electro-osmotic  drag  coefficient.  h2o  is 
the  diffusion  coefficient  of  water  in  the  membrane,  Am?p  is  the 
permeability  and  /xH2qi  the  viscosity  of  liquid  water.  Hence,  the 
transport  of  water  through  the  membrane  is  given  by  electro- 
osmotic  drag  and  diffusion  due  to  a  concentration  and  pressure 
gradient  across  the  membrane. 


AqLQff^~^Cjdt(Cqj  Tq) 
i 

=  -EC^%TeCfifout-<r^)+^effA?,sg^sg(7;-7i) 

i 

— A^vapiVg,  phase-  (17) 

Rearranging  Eq.  (17)  results  in 

-FcfP 

i 

_  1  r  /  ArOUt,rec Trec  _  at in,rec Trec  \ 

A  4>°ut  nq,i  1  q, ini 

+^t4g_  _-+ 

A q  Aq 

Integrating  Eq.  (7)  along  the  channel  length  yields 

rLe  ff  rLe  ff 


j 


dt[cc,o2(x’  0]  dx 


-/ 


dANc.o  Jx,  t)]  dx 


dywq  rL* 

■~i, FJo  1{x-,)ix- 


Using  the  mean  value  theorem  integration  leads  to 


(IB) 


(19) 


4.  Discretization 


The  stack  model  presented  here  simulates  the  most  important 
parameters  for  the  operation  of  a  PEMFC  stack,  which  are  illus¬ 
trated  in  Fig.  1.  To  derive  appropriate  and  rapidly  computable 
equations,  the  system  of  partial  differential  equations  for  the 
REVs  described  in  Section  3  is  integrated  along  the  effective 
channel  length  Le ff.  This  corresponds  to  the  second  step  of  the 
model  derivation  as  described  in  Section  2.  The  model  equations 
after  step  one  of  the  model  derivation  are  formulated  along  the 
channel,  that  is,  in  the  x-direction.  This  corresponds  to  assum¬ 
ing  straight  channels.  Thus,  Le ff  is  the  geometrical  length  of 
the  channel  measured  along  the  channel,  even  in  a  serpentine 
flow-field.  With  the  aid  of  the  Gaussian  law  and  the  mean  value 
theorem,  ordinary  differential  equations  (ODEs)  in  time  were 
derived  which  are  discretized  to  the  inlet  and  outlet  of  the  gas 
channels  or  the  cell’s  solid  material,  respectively.  In  the  follow¬ 
ing  the  approach  is  exemplified  for  Eqs.  (5)  and  (7).  Integration 
of  Eq.  (5)  along  the  channel  length  yields 


E/^eff 

CiAqdt  /  (Cqq(x,t)Tq(x,t))dx 
i  70 

ErLc  ff 

Ci  /  dx(Nq  i(x,  t)Tq(x,  t))  dx 

i  7o 

rLe  ff 

3~Aq  SgUSg  /  (Ts(x,  t) 

Jo 

ff 

1  /  phase  (-U 

Jo 


-Tq(x,  t))  dx  —  A Hy 


t)  dx. 


(16) 


dtcc,o2(t) 


1 


^cT'eff 


[AfXW  -  <ofC(0] 


dyWgjjt) 

AC4F 


(20) 


Thus,  balance  equations  are  derived  for  the  REC  which  are 
reduced  to  net  flows  at  the  channel  inlet  and  outlet.  Under  the 
assumption  that  all  gas  channels  in  a  cell  behave  like  the  REC 
net  flows  at  the  cell  inlet  and  outlet  are  calculated  through  mul¬ 
tiplication  of  the  channel  net  flows  with  the  number  of  channels 
in  each  electrode.  Integration  of  Eq.  (6)  directly  yields  a  net 
energy  balance  equation  for  the  solid  material  of  one  complete 
cell  in  the  stack.  Hence,  a  system  of  equations  is  derived,  which 
describes  a  complete  cell  module. 

In  step  three  of  the  model  derivation  the  stack  is  modeled  as 
several  coupled  cell  modules.  To  preserve  computational  effi¬ 
ciency  we  assume  that  each  cell  in  the  stack  works  identically 
and  that  the  fuel  is  distributed  equally  among  the  cells.  Fig.  4 
gives  an  overview  of  the  input  and  solution  variables  of  the  model 
and  illustrates  how  the  variables  are  adapted  from  cell  to  stack 
level. 


5.  Numerical  solution  method 

The  discretized  model  equations  form  a  linear-implicit  sys¬ 
tem  of  differential-algebraic-equations  (DAEs),  which  can  be 
written  as 

M(y,  t)y  =  f(y,  t)  (21) 

with  a  singular  mass  matrix  M(y,  t).  The  system  is  of  index  1, 
that  is,  one  derivation  is  necessary  to  transform  the  DAE-system 
into  a  system  of  ODEs.  Numerically  solving  a  DAE-system  is 
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Input  variables  Solution  variables 


i=H2,  h2o\  h2o' 

j=02,H20\H20' 


Fig.  4.  Overview  of  the  boundary  conditions  and  the  solution  variables  of  the 
model.  The  conversion  of  the  extensive  variables  from  the  stack  scale  to  the  cell 
scale  is  done  by  division  through  the  number  of  cells  nceiis  or  by  multiplication 
with  n  Ceiis,  respectively. 

generally  more  complex  than  solving  ODEs,  as  DAE-systems 
require  consistent  initial  conditions.  The  model  was  imple¬ 
mented  in  MATLAB™  using  the  built-in  solver  ode  15s  of  the 
ODE  Solver  environment  which  demonstrated  good  stability. 
Ode  15s  is  a  variable-order  solver  based  on  the  numerical  dif¬ 
ferentiation  formulas  capable  of  solving  stiff  DAE-systems  of 
index  1  [40]. 

6.  Parameter  identification 

The  mathematical  model  described  contains  a  high  num¬ 
ber  of  constants  and  parameters,  which  can  be  grouped 
into  stack-dependent  geometrical,  electrochemical  and  physi¬ 
cal  parameters,  as  well  as  stack-independent  physical  constants. 
Finding  suitable  values  for  these  parameters  is  one  of  the  crucial 
points  for  the  success  of  the  model.  The  identification  methods 
exemplified  for  the  stack  shown  in  Fig.  5  are  easily  transferable 
to  different  stack  designs. 

Most  of  the  geometrical  stack  and  flow-field  parameters  can 
be  taken  from  direct  measurements  as  illustrated  in  Fig.  2.  How¬ 
ever,  due  to  the  simplified  model  structure  several  geometrical 


Fig.  5.  Fuel  cell  stack  used  for  parameter  identification  and  model  validation. 


parameters  need  to  be  adjusted  accordingly.  The  stack  used 
for  model  validation  has  a  serpentine  flow-field.  As  the  model 
assumes  straight  channels  the  effective  channel  length  Leff  cor¬ 
responds  to  the  geometrical  channel  length  measured  along  the 
channel.  The  heat  exchange  area  between  solid  and  gas  A^sg 
equals  the  surface  of  one  channel  in  flow-field  q  standardized 
to  one  unit  length.  The  heat  exchange  area  between  solid  and 
surroundings  Ass  takes  into  account  the  surface  of  the  cells,  the 
endplates  and  the  cooling  fins.  The  surface  area  is  divided  by  the 
number  of  cells  and  standardized  to  one  unit  length.  The  scaling 
coefficient  dy  takes  into  account  the  enlargement  of  the  contact 
area  between  gas  and  membrane  due  to  diffusion  under  the  land 
areas  of  the  flow-field.  For  a  well-designed  flow-field  it  can  be 
assumed  that  the  whole  active  area  is  supplied  with  gases.  Thus, 
dy  equals  the  ratio  of  active  area  and  direct  contact  area  between 
gas  channels  and  GDL. 

The  physical  and  reaction  parameters  are  determined  by  mea¬ 
surements  and  literature  research  as  indicated  in  Table  1.  The 
cathodic  exchange  current  density  j'o?c,  the  symmetry  factor  a 
and  the  losses  in  open-circuit  voltage  ijoc  are  identified  through 
least-square-fits  of  Eq.  (1 1)  to  IV-curves.  The  heat  transfer  coef¬ 
ficient  between  solid  material  and  surroundings  Uss  is  identified 
on  the  basis  of  cooling  curves  monitoring  the  average  stack  tem¬ 
perature  rS)av  in  a  setting  without  gas  and  current  flow.  In  this 
case  the  following  equation  can  be  deduced  from  Eq.  (6): 

d^av  =  UssAs L(rsurW  _  T  (f))  (22) 

dt  psCsAs 
Integration  leads  to 

7i,av(0  =  (TSt av(fo)  -  TsM)  exp  ( +  7W.  (23) 

v  PsCsAs  J 

The  parameter  Uss  is  determined  through  least- square- fits  of  Eq. 
(23)  to  measured  cooling  curves. 

7.  Experiment 

In  order  to  identify  the  dynamic  behavior  of  a  portable  PEM 
fuel  cell  stack  an  experimental  investigation  was  carried  out. 
The  portable  stack  used  for  the  measurements  was  developed 
at  the  Fraunhofer  ISE.  It  consists  of  six  cells  with  an  active 
cell  area  of  30.2  cm2  each.  A  Gore  Primea  5510  MEA  with  a 
platinum  loading  of  0.4  mg  cm-2  is  used.  Cooling  fins  on  the 
outer  cells  allow  operation  with  passive  cooling  only.  Due  to 
its  small  geometrical  dimensions  the  stack  shows  a  pronounced 
dynamic  behavior  with  respect  to  stack  temperature  and  water 
management.  Therefore,  it  is  very  well  suited  for  a  validation 
study  of  the  dynamic  stack  model.  In  addition,  the  quick  response 
of  the  stack  to  changes  in  the  operating  conditions  allows  short, 
controlled  operation  under  critical  conditions  such  as  dry-out 
or  flooding  without  damaging  the  stack.  The  particular  stack  is 
used  for  the  implementation  of  NMPC  in  an  ongoing  project. 
Results  of  this  work  will  be  published  elsewhere. 

The  measurements  were  carried  out  at  a  computer-controlled 
test  stand,  which  controls  the  operating  conditions  of  the  stack  as 
well  as  the  data  acquisition.  Gas  flow  controllers  manage  the  flow 
of  hydrogen  into  the  anode  and  oxygen  or  air  into  the  cathode. 
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An  external  consumer  is  simulated  by  a  galvanostat,  which  can 
subject  the  stack  to  arbitrary  current  profiles.  The  gas  inlet  flows 
as  well  as  the  voltages  of  each  individual  cell  and  of  the  whole 
stack  are  measured.  The  stack  impedance  is  measured  at  1  kHz. 
The  stack  temperature  was  measured  with  four  thermocouples 
at  the  external  surface  of  the  stack  as  usually  done  in  integrated 
fuel  cell  systems.  As  indicated  in  Fig.  1  the  thermocouples  were 
placed  on  the  long  and  short  sides  of  the  top  and  of  one  of 
the  middle  cells  in  order  to  obtain  a  reasonable  average  value. 
It  is  assumed  that  due  to  the  high  thermal  conductivity  of  the 
bipolar  plates  the  temperature  distribution  throughout  the  stack 
balances  quickly.  This  assumption  is  supported  by  the  quick 
reaction  of  the  measured  stack  temperature  to  changes  in  the 
load.  In  addition,  a  comparison  of  the  measured  temperatures  of 
the  four  thermocouples  showed  a  maximum  deviation  of  only 
2K. 

The  stack  was  subjected  to  different  load  profiles,  e.g.  step 
and  jump  profiles  during  which  it  was  supplied  with  pure  hydro¬ 
gen  on  the  anode  side  and  dry  air  on  the  cathode  side.  For 
each  measurement  the  gas  flow  was  kept  at  a  constant  rate, 
which  corresponds  to  a  stoichiometry  of  2  for  the  highest  cur¬ 
rent  in  the  profile.  The  measurements  were  carried  out  under 
ambient  pressure.  The  temperature  of  the  surroundings  was 
Tsur  =  (298  zb  1)  K  throughout  the  experiment.  Each  profile  was 
measured  several  times  to  ensure  reproducibility. 

8.  Results  and  discussion 

For  the  validation  study  several  profiles  were  simulated.  Sta¬ 
ble  numerical  convergence  behavior  was  observed  in  all  cases. 
The  parameters  used  for  the  simulations  are  listed  in  Table  1 .  In 
the  following,  a  step  and  a  jump  profile  are  discussed.  For  the 
step  profile  stack  currents  of  I  =  1,  3,  6  A  were  applied.  The 
time- spans  for  each  step  were  3  min  on  the  ascent  and  6  min  on 
the  descent.  For  the  later  a  longer  interval  was  chosen  to  fur¬ 
ther  investigate  the  cooling  of  the  stack.  The  flow  rates  at  the 
gas  inlets  were  set  to  0.7  /  min-1  for  hydrogen  and  1.4 1  min-1 
for  air.  For  the  jump  profile  the  current  was  alternated  between 
1  A  and  5  A  with  a  hydrogen  flow  rate  of  0.6 1  min-1  and  an  air 
flow  rate  of  1 .2  /  min-1 .  The  simulations  were  carried  out  on  an 
AMD  Athlon  1533  MHz.  Short  computing  times  of  0.6  s  for  the 
step  profile  with  a  time  of  operation  of  23  min  and  0.8  s  for  the 
jump  profile  with  a  time  of  operation  of  30  min  are  proofs  of  the 
numerical  efficiency  of  the  model. 

8.1.  Model  validation 

Validation  of  the  model  is  performed  by  comparison  of  mea¬ 
sured  stack  temperatures  and  voltages  to  the  model  predictions. 
Figs.  6  and  7  show  a  comparison  between  measured  and  sim¬ 
ulated  stack  temperature  for  the  step  and  the  jump  profile, 
respectively.  The  measured  value  is  an  average  of  the  measure¬ 
ments  of  the  four  thermocouples.  Simulation  and  experiment 
show  good  correlation.  At  high  current  the  stack  temperature 
increases  due  to  the  different  loss  mechanisms  corresponding 
to  term  (b)  in  Eq.  (6).  Parametric  studies  showed  that  activa¬ 
tion  losses  at  the  cathode  are  the  dominant  loss  mechanism.  The 


Fig.  6.  Comparison  between  simulated  and  measured  stack  temperature  for  a 
current  step  profile.  Simulation  and  experiment  show  good  correlation.  The  break 
on  step  4  is  caused  by  the  phase  change  enthalpy  of  water. 


cooling  of  the  stack  is  dominated  by  heat  exchange  between 
solid  material  and  surroundings  (term  (d)  in  Eq.  (6)),  whereas 
only  a  small  amount  of  heat  is  dissipated  by  the  exhaust  gases. 
It  is  noticeable  that  current  changes  affect  the  stack  temper¬ 
ature  almost  without  time-delay,  even  though  the  temperature 
was  measured  at  the  outside  of  the  stack  as  indicated  in  Fig.  1. 
This  is  explained  by  the  small  thermal  mass  and  the  high  heat 
conductivity  of  the  solid  stack  material,  hence  validating  the 
assumption  of  a  homogeneous  stack  temperature  in  the  y-  and 
z-direction.  Although  the  simulated  and  measured  stack  tem¬ 
peratures  agree  well,  there  are  also  slight  deviations.  It  can  be 
noted,  that  the  deviations  are  mostly  due  to  a  slower  increase 
and  a  faster  decrease  of  the  simulated  stack  temperature  at  low 
current  (Fig.  6,  steps  1  and  2;  Fig.  7,  steps  with  odd  numbers), 
whereas  at  high  current  the  slopes  agree  well.  As  the  water  pro¬ 
duction  and  therefore  the  degree  of  humidification  of  membrane 
and  electrode  is  correlated  with  the  current  density,  it  can  be 
concluded  that  a  humidification-dependent  heat  source  is  not 
modeled  comprehensively.  This  can  be  the  heat  dissipation  due 
to  insufficient  humidification  of  the  membrane-electrode  assem¬ 
bly  as  well  as  inhomogeneous  current  distributions  and  losses 
throughout  the  stack.  The  break  on  step  4  in  Fig.  6  is  caused 


Time  [min] 

Fig.  7.  Comparison  between  simulated  and  measured  stack  temperature  for  a 
current  jump  profile.  Simulation  and  experiment  agree  well.  Deviations  can  be 
explained  by  heat  dissipation  due  to  insufficient  humidification  of  the  electrode 
and  inhomogeneous  current  density  distribution. 
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Fig.  8.  Comparison  between  simulated  and  measured  stack  voltage  for  a  current 
step  profile.  A  good  correlation  is  obtained.  Deviations  between  model  predic¬ 
tions  and  experiment  are  connected  with  critical  operating  conditions  of  the 
stack. 
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Fig.  10.  Simulation  results  for  the  step  profile  indicate  a  high  concentration  of 
liquid  water  on  the  cathode  side  for  steps  2-4.  This  coincides  with  low  values  of 
simulated  and  measured  impedance  on  these  steps.  Observation  of  the  average 
concentration  of  liquid  water  and  the  stack  impedance  can  be  used  to  predict 
critical  operating  conditions  like  flooding  or  dry-out. 


by  phase  change  enthalpy  of  water.  In  the  simulation  all  liquid 
water  has  evaporated  at  this  point  as  shown  in  Fig.  1 1 .  Heat  con¬ 
sumption  due  to  evaporation  stops,  which  leads  to  an  increase  of 
the  stack  temperature.  In  the  model  this  is  reflected  by  a  change 
of  sign  of  term  (a)  in  Eq.  (6)  (A^, phase  is  negative,  when  evapo¬ 
ration  takes  place  and  positive  elsewise).  In  reality  a  remainder 
of  water  will  longer  exist  in  the  GDL  or  the  membrane  leading 
to  a  milder  change-over  between  the  two  regimes.  A  detailed 
discussion  on  water  balance  is  given  below.  Similar  breaks  indi¬ 
cating  the  switch  between  the  two  regimes  are  visible  on  steps 
5,  7  and  9  in  Fig.  7.  As  the  temperature  level  increases  with  the 
step  number  the  amount  of  accumulated  water  throughout  the 
preceding  high  current  step  gets  smaller  and  evaporation  gets 
faster.  Therefore,  the  breaks  on  higher  numbered  steps  occur 
closer  to  the  preceding  load  change.  As  the  stack  was  dry  at 
the  beginning  of  the  profile  no  breaks  occur  on  steps  1  and  3. 
The  validity  of  this  interpretation  is  shown  after  the  discussion 
of  water  balance  below.  A  comparison  between  measured  and 
simulated  stack  voltage  is  shown  in  Figs.  8  and  9.  Simulated  and 
measured  stack  voltage  agree  well.  Deviations  between  model 
predictions  and  experiment  are  connected  with  the  water  bal¬ 
ance  of  the  stack  as  discussed  below  and  in  Fig.  12.  A  satisfying 
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Fig.  9.  Comparison  between  simulated  and  measured  stack  voltage  for  a  current 
jump  profile.  Measurement  and  experiment  agree  well.  The  model  parameters 
remain  unchanged  for  each  case  of  the  validation  study. 


convergence  behavior  was  also  achieved  for  the  other  load  pro¬ 
files  in  the  validation  study.  Thus,  the  model  predictions  for  the 
stack  temperature  and  the  voltage  are  validated. 

8.2.  Water  balance  and  critical  operating  conditions 

Water  management  plays  an  important  role  for  the  opera¬ 
tion  of  a  PEMFC.  Low  concentrations  of  liquid  water  lead  to  a 
decrease  in  membrane  conductivity  whereas  an  accumulation  of 
liquid  water  can  block  the  gas  flow  due  to  flooding  of  electrode 
and  gas  channels.  The  stack  model  presented  here  incorporates 
mass  balances  of  water  vapor  and  liquid  water  (Eqs.  (8)  and  (9)), 
which  are  linked  through  a  phase  change  term  (Eq.  (10)).  To 
allow  prediction  of  critical  operating  conditions  in  an  integrated 
fuel  cell  system,  the  stack  model  predicts  the  average  concen¬ 
trations  of  liquid  water  at  the  anode  and  cathode  side  and  the 
impedance  of  the  stack.  Fig.  10  shows  corresponding  simulation 
results  and  experimental  values  of  the  stack  impedance  for  the 
step  profile.  Simulation  results  indicate  a  high  concentration  of 
liquid  water  at  the  cathode  side  for  steps  2-4.  This  coincides  with 
low  values  of  simulated  and  measured  impedance  on  steps  2-4. 
As  the  membrane  model  is  steady-state,  the  simulated  values 
change  rapidly  upon  a  load  change,  whereas  the  experimen¬ 
tal  results  indicate  slower  variations  of  membrane  conductivity. 
However,  the  simulated  impedance  reflects  the  trend  of  the  mea¬ 
surement  well.  In  order  to  cover  the  dynamic  behavior  of  the 
water  impregnated  into  the  membrane,  a  dynamic  membrane 
model  is  necessary.  Yet,  models  that  describe  the  dynamic  water 
transport  through  the  membrane  accurately  are  usually  com¬ 
putationally  expensive,  see  e.g.  Ref.  [35].  The  development  of 
a  computationally  efficient,  dynamic  membrane  model  is  very 
challenging,  but  is  seen  as  one  of  the  key  points  to  improve  the 
predictive  capability  of  dynamic  PEMFC  models.  To  illustrate 
the  evolution  of  the  average  concentration  of  liquid  water  in  the 
cathode,  Fig.  1 1  shows  the  relation  between  phase  transition  of 
water,  its  partial  pressure  and  the  saturation  pressure  for  the  step 
profile.  Water  condenses  when  pu2ow  >  ^c,  sat  and  liquid  water 
evaporates  when  pu2ow  <  ^c,  sat  (Eq.  (10)).  It  is  assumed  that 
the  stack  does  not  contain  liquid  water  at  the  beginning  of  the 
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Fig.  11.  Simulation  results  visualizing  the  relation  between  phase  transition  of 
water,  its  partial  pressure  pc(H2Ov)  and  the  saturation  pressure  PCj  sat  of  water 
vapor  on  the  cathode  for  the  step  profile.  Condensation  occurs  between  point  a 
and  b  as  pc( H20v)  is  higher  than  Pc?  sat.  Liquid  water  evaporates  between  points 
b  and  c  as  pc(H2Ow)  <  PC;  sat  holds. 


measurement.  Condensation  occurs  during  steps  2  and  3  starting 
at  point  (a),  where  pu2ow  is  higher  than  Pc?  sat  due  to  the  high 
production  of  water  in  the  electrochemical  reaction.  This  coin¬ 
cides  with  an  increasing  concentration  of  liquid  water  in  these 
steps  in  Fig.  10.  Condensation  ends  as  pu2ow  decreases  due  to 
lower  water  production  after  the  load  reduction  (b).  Existing  liq¬ 
uid  water  evaporates  in  the  course  of  several  minutes  (c).  After 
this  point  due  to  the  higher  temperatures  and  the  low  load  water 
exists  in  vapor  form  only  and  condensation  or  evaporation  does 
not  take  place. 

Predicting  the  impedance  and  the  average  concentration  of 
liquid  water,  the  model  allows  the  early  detection  of  critical  oper¬ 
ating  conditions  such  as  flooding  or  dry-out.  However,  the  model 
does  not  yet  predict  the  corresponding  quantitative  decrease  in 
the  stack  voltage  as  the  two-phase  effects  require  detailed  mod¬ 
els  on  the  micro-scale,  which  are  not  agreeable  to  the  demand 
of  limited  computational  effort.  Fig.  12  illustrates  the  effects  of 
the  water  balance  on  the  stack  voltage.  The  measured  voltage 
curve  shows  three  characteristics:  humidification,  flooding  and 
dry-out.  At  step  2  the  stack  voltage  increases  as  the  membrane 


Fig.  12.  Comparison  of  measured  and  simulated  voltage  for  the  current  step 
profile  indicating  the  effects  of  the  water  balance.  The  measured  voltage  curve 
shows  three  characteristics:  humidification,  flooding  and  dry-out.  The  simulated 
concentration  of  liquid  water  in  Fig.  10  indicates  the  breakdown  of  the  stack 
voltage  due  to  flooding. 


and  electrode  take  up  water.  The  opposite  effect  is  observed  dur¬ 
ing  steps  4  and  5.  The  simulated  average  concentration  of  liquid 
water  in  Fig.  10  as  well  as  the  impedance  indicate  a  low  concen¬ 
tration  of  water  in  the  stack.  In  step  3  the  stack  voltage  shows 
flooding  effects,  which  are  not  captured  by  the  model.  How¬ 
ever,  the  simulated  concentration  of  liquid  water  indicates  the 
breakdown  of  the  stack  voltage  due  to  flooding.  Thus,  observ¬ 
ing  the  characteristics  of  the  voltage  allows  a  conclusion  on  the 
water  balance  of  the  stack.  This  also  supports  the  validity  of 
the  interpretation  of  Fig.  7.  During  step  1  in  Fig.  9  the  voltage 
increases  rapidly  indicating  a  low,  but  increasing  degree  of  mem¬ 
brane  and  electrode  humidification,  whereas  throughout  steps  3, 
5,  7  and  9  the  voltage  decreases  due  to  drying-out  of  the  mem¬ 
brane.  For  higher  numbered  steps  the  decrease  becomes  much 
faster.  Although  the  model  does  not  capture  the  voltage  decrease 
quantitatively,  the  position  of  the  breaks  in  the  simulated  stack 
temperature  in  Fig.  7  shows  that  the  regimes  of  high  and  low 
water  concentration  are  modeled  appropriately. 

As  Fig.  12  indicates  the  predictive  capability  of  the  model 
shows  limitations  when  the  stack  is  operated  in  a  high  cur¬ 
rent  range.  This  is  due  to  the  negligence  of  the  GDL  and  the 
simple  model  approach  used  for  the  membrane.  Therefore,  the 
highly  complex  two-phase  effects  in  the  membrane  and  the  GDL 
are  not  covered.  To  develop  a  computationally  efficient  model, 
that  includes  the  corresponding  phenomena  goes  far  beyond  the 
scope  of  the  present  paper.  In  addition,  the  predictive  capability 
of  the  model  might  show  limitations  for  PEMFC  stacks  with  a 
strong  spatial  inhomogeneity  of  temperature  or  electric  poten¬ 
tial.  This  is  due  to  the  fact  that  the  gradients  in  temperature  and 
electric  potential  across  the  stack  need  to  be  small  to  ensure  the 
validity  of  the  model  approach. 

9.  Conclusions 

In  this  work,  a  numerically  efficient  model  was  developed  to 
study  the  dynamic  response  of  a  PEMFC  stack  subjected  to  load 
changes.  The  attribute  of  numerical  efficiency  was  imposed  in 
order  for  the  model  to  be  applicable  for  NMPC.  Therefore,  a  rea¬ 
sonable  compromise  was  found  between  computational  effort 
and  physical  accuracy.  Mainly  three  model  properties  make 
our  calculations  computationally  efficient:  firstly,  only  one  cell 
is  simulated  and  the  results  are  then  upscaled  to  account  for 
the  whole  stack.  Secondly,  the  spatial  dependency  of  the  fuel 
cell  parameters  is  neglected  as  the  energy  and  mass  balance 
equations  on  which  the  model  is  based  are  reduced  to  ordi¬ 
nary  differential  equations  in  time.  And  thirdly,  the  GDL  is  not 
modeled  and  the  membrane  is  assumed  to  be  at  steady-state  as 
a  detailed  modeling  of  these  cell  parts  would  require  massive 
computational  effort.  As  a  result  of  this  approach  the  model 
proves  to  be  computationally  efficient:  an  operating  time  of  0.5 
h  can  be  simulated  in  less  than  1  s.  Despite  the  relatively  simple 
model  approach,  the  model  and  the  calculation  method  could  be 
validated  against  experimental  data  for  different  realistic  load 
profiles  through  a  comparison  of  measured  and  simulated  stack 
voltage  and  stack  temperature  on  time  scales  of  up  to  30  min. 
The  model  was  used  to  study  the  run  of  the  stack  voltage  and  the 
stack  temperature  curves.  Parametric  studies  showed  that  activa- 
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tion  losses  at  the  cathode  cause  the  main  mechanism  of  heating, 
whereas  the  heat  exchange  between  solid  and  surroundings  is 
the  main  mechanism  for  cooling.  Slight  deviations  of  the  slope 
of  the  simulated  stack  temperature  at  low  current  indicate  that 
a  humidification-dependent  heat  source  is  not  modeled  com¬ 
prehensively,  which  can  be  heat  dissipation  due  to  insufficient 
humidification  of  the  membrane-electrode  assembly  as  well  as 
inhomogeneous  current  distributions  throughout  the  stack.  Fur¬ 
thermore,  it  is  mentionable  that  current  changes  affect  the  stack 
temperature  almost  without  time-delay,  although  the  tempera¬ 
ture  was  measured  at  the  outside  of  the  stack.  This  is  explained 
by  the  small  thermal  mass  and  the  high  heat  conductivity  of  the 
solid  material. 

In  order  to  ensure  reliable  operation  of  a  PEMFC  under  load 
changes  it  is  essential  to  prevent  critical  operating  conditions  of 
flooding  and  dry-out,  which  depend  on  the  highly  complex  phe¬ 
nomena  of  water  balance  of  the  stack.  Therefore,  we  put  special 
emphasis  on  the  study  of  possibilities  and  limitations  of  the  water 
balance  description  with  respect  to  critical  operating  conditions 
within  such  a  reduced  model.  Our  model  incorporates  mass  bal¬ 
ance  equations  for  water  vapor  and  liquid  water  including  a  term 
for  the  phase  change.  For  the  purpose  of  studying  the  capabil¬ 
ity  of  this  approach  the  stack  was  operated  in  critical  states  of 
flooding  and  dry-out.  Through  a  comparison  of  measured  and 
simulated  stack  voltage  three  characteristics  of  the  stack  voltage 
indicating  humidification,  flooding  and  dry-out  of  the  stack  are 
identified.  Thereby  indicating  how  the  observation  of  the  stack 
voltage  allows  a  conclusion  on  the  water  balance  of  the  stack. 
The  model  was  successfully  applied  to  predict  the  regimes  of 
high  and  low  water  concentrations  in  the  stack.  It  is  shown,  that 
by  monitoring  the  simulated  stack  impedance  and  the  average 
concentration  of  liquid  water  in  the  stack  critical  states  of  oper¬ 
ation  can  be  predicted.  However,  the  model  does  not  yet  capture 
the  corresponding  quantitative  change  of  the  stack  voltage  as 
the  underlying  two-phase  effects  require  detailed  descriptions 
on  the  micro-scale  which  are  not  compatible  with  the  demand 
for  limited  computing  time. 

With  our  model  we  have  demonstrated  that  the  dynamic 
response  of  a  PEMFC  stack  under  load  changes  can  be  sim¬ 
ulated  in  short  computing  times  with  sufficient  accuracy.  Efforts 
are  currently  underway  to  implement  NMPC  algorithms  based 
on  this  model.  In  order  to  further  improve  the  predictive  capabil¬ 
ities  of  the  model  future  work  could  include  the  development  of 
a  computationally  efficient,  dynamic  membrane  model  as  well 
as  the  integration  of  a  GDL  model  to  improve  simulation  results 
when  the  stack  is  operated  in  a  high  current  range. 
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Appendix  A 


A.l.  Supporting  equations 


The  internal  energy  Ukj  of  species  i  in  volume  element  k  is 
given  by 


UkJ  =  CkjCiTk. 


(24) 


The  open-circuit  voltage  Uoc  is  calculated  using  the  Nernst- 
equation: 

Uoc  =  u°0c  +  ^-(T  -  T° )  -  In  (^)  ,  (25) 

where  U°c  is  the  open-circuit  voltage  at  reference  conditions, 
ASmoi  is  the  molar  entropy  of  the  overall  reaction  and  T°  is  the 
reference  temperature.  V/  denotes  the  stoichiometry  factor,  pi  is 
the  partial  pressure  and  p ?  is  the  partial  pressure  of  species  i  at 
reference  conditions.  The  partial  pressure  pi  of  the  gas  species 
i  is  calculated  using  Dalton’s  law: 

Pi  =  ^~Pq,  (26) 

/  ;N<lJ 
j 


where  Pq  is  the  average  pressure  in  gas  channel  q.  Assuming  that 
the  gases  in  the  channels  can  be  treated  as  ideal,  the  following 
equation  is  used  to  calculate  the  concentration  cq ?  /  from  the  molar 
flow  Nqj : 


c 


q,i  — 


Nq,j 

AqVqi 


(27) 


where  vqj  is  the  velocity  of  species  i  in  electrode  q. 

The  electro-osmotic  drag  coefficient  n^g  is  calculated 
according  to  Springer  et  al.  [34]: 


ttdrag  =  (5/44)V.  (28) 

The  following  empirical  equation  is  used  to  calculate  the  satu¬ 
ration  pressure  of  water  vapor  [34] : 

Pq  sat=c pixl()(cP2+cp3(Tq-213)+cp4(Tq-213)2+cp5(Tq-273)3)  (29) 

with  cp\  =  101,  325,  cp2  =  —2.18,  c/>3  =  2.95  x  10-2,  cp4 
=  —9.18  x  10~5  andcps  =  1.44  x  10-7.  The  diffusion  coef¬ 
ficient  of  water  in  the  membrane  Dm? h20  as  a  function  of  A,m 
and  Ts  was  described  by  Golbert  and  Lewin  [36]  as 

An,H20  =  «drag£>^fH20  exp  \cD\  ^  -d—  -  2^  J  (30) 

with  the  empirical  values  cd i  =  2416  and  7^, ref  =  303  K  [36]. 
^mH2o  diffusion  coefficient  of  water  in  the  membrane  at 
reference  conditions. 

The  activity  of  water  in  channel  q  is  calculated  according  to 
Golbert  and  Lewin  [36]: 


<2h2ov 


Nq,  H2Ov  Pq 

^2Nq,i  ^?’sat 


(31) 
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The  reaction  entropy  of  the  overall  reaction  AS^i  for  the  con¬ 
sumption  of  1  mol  H2  is  calculated  from  the  reaction  entropy  at 
the  anode  A Sa  and  the  cathode  A Sc: 

A  Sc 

ASmol  =  ASa+^.  (32) 
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